function mu_new=mu(mu_old,t,f)
if(abs(mu_old(1,3))<0.99999)
    mu_new(1,1)=(sin(t)/sqrt(1-mu_old(1,3)^2))*(mu_old(1,1)*mu_old(1,3)*cos(f)-mu_old(1,2)*sin(f))+mu_old(1,1)*cos(t);
    mu_new(1,2)=(sin(t)/sqrt(1-mu_old(1,3)^2))*(mu_old(1,2)*mu_old(1,3)*cos(f)+mu_old(1,1)*sin(f))+mu_old(1,2)*cos(t);
    mu_new(1,3)=-sin(t)*cos(f)*sqrt(1-mu_old(1,3)^2)+mu_old(1,3)*cos(t);
else(mu_old(1,3)>0.99999)
    mu_new(1,1)=sin(t)*cos(f);
    mu_new(1,2)=sin(t)*sin(f);
    mu_new(1,3)=mu_old(1,3)/abs(mu_old(1,3))*cos(f);
end
    